% Copyright (C) 2009,2010,2011,2012  Marco Restelli
%
% This file is part of:
%   LDGH -- Local Hybridizable Discontinuous Galerkin toolkit
%
% LDGH is free software: you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% LDGH is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
% or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
% License for more details.
%
% You should have received a copy of the GNU General Public License
% along with LDGH. If not, see <http://www.gnu.org/licenses/>.
%
% author: Marco Restelli                   <marco.restelli@gmail.com>


function uuI = dgcomp_fast_output_intrp(grid,base,data,x)
% uuI = dgcomp_fast_output_intrp(grid,base,data,x)
%
% Interpolate a time series in an arbitrary point (provided it is
% contained in an element for which the time series is recorded).
%
% grid: the computational grid
% base: FE basis
% data: the high frequency data structure returned by
%       collect_fast_output
% x: the point where the interpolation is required (column array)

 d = grid.d;

 % 1) locate the correct element
 toll = 1e-10;
 found = 0;
 ie = 1;
 while( and( not(found) , ie<=length(data.ie) ) )
   e = grid.e(data.ie(ie));
   xi(2:d+1) = e.bi*(x-e.x0); % barycentric coords
   xi(1) = 1 - sum(xi(2:d+1));
   ximin = min(xi);
   if(ximin>=-toll)
     found = 1;
   else
     ie = ie+1;
   end
 end

 if(not(found))
   error('Can not locate the requested point in any of the time history elements')
 end

 % 2) evaluate the basis functions
 uPHI = zeros(ubase.pk,1);
 for i=1:ubase.pk
   uPHI(i) = ev_pol(ubase.p_s{i},xi(2:end)');
 end
 pPHI = zeros(pbase.pk,1);
 for i=1:pbase.pk
   pPHI(i) = ev_pol(pbase.p_s{i},xi(2:end)');
 end

 % 3) interpolate the time history
 uup = zeros(d+1,length(data.times));
 for i=1:d
   for j=1:ubase.pk
     uup(i,:) = uup(i,:) + ...
       permute(data.uu(i,j,ie,:),[1 4 2 3])*uPHI(j);
   end
 end
 for j=1:pbase.pk
   uup(d+1,:) = uup(d+1,:) + ...
       permute(data.pp(j,ie,:),[1 3 2])*pPHI(j);
 end

return
